{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Catching close encounters using exceptions\n",
    "Sometimes one is interested in catching a close encounter between two planets. This can easily be done with REBOUND. What you do when a close encounter happens is up to you.\n",
    "\n",
    "Some integrators are better suited to simulate close encounters than others. For example, the non-symplectic integrator IAS15 has an adaptive timestep scheme that resolves close encounters very well. Integrators that use a fixed timestep like WHFast are more likely to miss close encounters.\n",
    "\n",
    "Let's start by setting up a two-planet system that will go unstable on a short timescale:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rebound\n",
    "import numpy as np\n",
    "def setupSimulation():\n",
    "    sim = rebound.Simulation()\n",
    "    sim.integrator = \"ias15\" # IAS15 is the default integrator, so we don't need this line\n",
    "    sim.add(m=1.)\n",
    "    sim.add(m=1e-3,a=1.)\n",
    "    sim.add(m=5e-3,a=1.25)\n",
    "    sim.move_to_com()\n",
    "    return sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's integrate this system for 100 orbital periods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = setupSimulation()\n",
    "sim.integrate(100.*2.*np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rebound exits the integration routine normally. We can now explore the final particle orbits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<rebound.Orbit instance, a=4.725537568103762 e=0.7145243085874022 inc=0.0 Omega=0.0 omega=2.7372774231509887 f=4.39210305882742>\n",
      "<rebound.Orbit instance, a=1.0429076415181782 e=0.12110823999615411 inc=0.0 Omega=0.0 omega=5.927708717019177 f=4.786502097969678>\n"
     ]
    }
   ],
   "source": [
    "for o in sim.orbits():\n",
    "    print(o)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the orbits of both planets changed significantly and we can already speculate that there was a close encounter.\n",
    "\n",
    "Let's redo the simulation, but this time set the `sim.exit_min_distance` flag for the simulation. If this flag is set, then REBOUND calculates the minimum distance between all particle pairs each timestep. If the distance is less than `sim.exit_min_distance`, then the integration is stopped and an exception thrown. Here, we'll use the [Hill radius](https://en.wikipedia.org/wiki/Hill_sphere) as the criteria for a close encounter. It is given by $r_{\\rm Hill} \\approx a \\sqrt{\\frac{m}{3M}}$, which is approximately 0.15 AU in our case. \n",
    "\n",
    "This setup allows us to catch the exception and deal with it in a customized way.  As a first example, let's catch the exception with a `try`-`except` block, and simply print out the error message. Additionally, let's store the particles' separations while we're integrating:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Two particles had a close encounter (d<exit_min_distance).\n"
     ]
    }
   ],
   "source": [
    "sim = setupSimulation() # Resets everything\n",
    "sim.exit_min_distance = 0.15\n",
    "Noutputs = 1000\n",
    "times = np.linspace(0,100.*2.*np.pi,Noutputs)\n",
    "distances = np.zeros(Noutputs)\n",
    "ps = sim.particles # ps is now an array of pointers. It will update as the simulation runs.\n",
    "try:\n",
    "    for i,time in enumerate(times):\n",
    "        sim.integrate(time)\n",
    "        dp = ps[1] - ps[2]   # Calculates the coponentwise difference between particles \n",
    "        distances[i] = np.sqrt(dp.x*dp.x+dp.y*dp.y+dp.z*dp.z)\n",
    "except rebound.Encounter as error:\n",
    "    print(error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Encounter` does currently not tell you which particles had a close encounter. But you can easily search for the pair yourself (see below). \n",
    "\n",
    "Here, we already know which bodies had a close encounter (the two planets) and we can plot their separation as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAAHACAYAAACVhTgAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB5SElEQVR4nO3dd3hUZdoG8PvMTCa990oSWggphB5AioCAiICoKCoI2LGi68q6a1tX3V1dy6pgo1hBpShFlN47BBJICJCEhJAe0pOp5/tjkqx8AilM8k65f9c1l2Qy5Q6GmXnO+5znlWRZlkFERERERERXpRAdgIiIiIiIyNKxcCIiIiIiImoBCyciIiIiIqIWsHAiIiIiIiJqAQsnIiIiIiKiFrBwIiIiIiIiagELJyIiIiIiohawcCIiIiIiImqBSnSAzmY0GnHx4kW4u7tDkiTRcYiIiIiISBBZllFdXY2QkBAoFNdeU7K7wunixYsIDw8XHYOIiIiIiCxEXl4ewsLCrnkbuyuc3N3dAZj+cjw8PASnISIiIiIiUaqqqhAeHt5cI1yL3RVOTe15Hh4eLJyIiIiIiKhVp/BwOAQREREREVELWDgRERERERG1gIUTERERERFRC1g4ERERERERtYCFExERERERUQtYOBEREREREbWAhRMREREREVELWDgRERERERG1gIUTERERERFRC1g4ERERERERtYCFExERERERUQtYOBEREREREbWAhRMREREREVELVKIDEJF5nS6sxroTF7E/qwyRvq4Y1t0Pw7r5wdfNUXQ0IiIiIqvFwonIBpwrqcG64wVYd+IizhTXNF9/KOcSfjhyAQAQG+yBG7r7YVh3PwyI9IGTg1JUXCIiIiKrI8myLIsO0Zmqqqrg6emJyspKeHh4iI5D1G555XVYe+Ii1h0vwKmCqubrHZQSRvTwx40xgcgurcGuM6XIKKy+7L5qlQIDIr0xrJs/bujuh9hgDygUUmf/CERERERCtaU2YOFEZEUKKuux/kQB1p4owPG8iubrlQoJQ7v5YVJCMG7qHQRPZ4fL7ldSrcGes6XYfbYUu8+UorCq4bLv+7iqMaSrL4Z1M61IhXm7dMaPQ0RERCQUC6drYOFE1qa4ugG/pBZi3YmLOJRzqfl6hQQMjvbFLQkhGB8XBB9XdaseT5ZlnCsxrUTtPlOK/VllqNUaLrtNlJ8rhnXzww3d/XBjTABUSs6RISIiItvDwukaWDiRNSiv1WJjWmHzkAfj7/6VDoj0xi0JIZgQH4QAd6frfi6dwYiUvIrGQqoExy9UwvC7J0yO9sXCe/vCy6V1hRkRERGRtWDhdA0snMjSrTtxEfO/Pw6t3th8XWK4FyYlBOPm+GCEeDl36PNXNeiw71wZ9pwtxcojF1CrNSDKzxVfzOqPaH+3Dn1uIiIios7EwukaWDiRJduYVoB53x6DwSgjJsgdk/uE4paEYIT7iDnnKKOwCnOXHkZ+RT08nR2w8N6+GNLVT0gWIiIiInNrS23AExeILMSmU0V4vLFouq1vKDY8eQMeHdlVWNEEADFBHlgzbyiSIrxQWa/DzC8OYvnBXGF5iIiIiERh4URkAbZlFOOxb45Ab5QxuU8I/n17osWMB/d3d8R3Dw7GrYkh0BtlvLAqFf9Yf+qy86CIiIiIbB0LJyLBdmaW4OGvj0BnkDExPhjv3JEIpYUUTU2cHJR4/64+eGZMDwDAZ7uy8fBXR1Cr0QtORkRERNQ5WDgRCbT3bCke/PIwtHojbooNxHt39bHY0d+SJOGpMd3xwd1JUKsU2JxehNsX7cPFinrR0YiIiIg6nGV+QiOyAweyyjB32WFo9EaMjgnAhzP6wsFCi6bfuzUxBMsfGgw/NzXSC6ow+aM9SPndZrxEREREtsjyP6UR2aAj58sxe+kh1OsMGNHDHx/f2xdqlfX8c+wb4Y0184YiJsgdJdUaTP9kH9afKBAdi4iIiKjDWM8nNSIbkZJXgVmLD6FOa8Cwbn745L5+cFQpRcdqszBvF/z46BDcGBMAjd6Ied8exX+3nIGd7XBAREREdoKFE1EnSr1Qifu+OIAajR6Do33w2cz+cHKwvqKpiZujCp/N7I+5w6IAAO9sysT8749DozcITkZERERkXiyciDrJqYtVuPeLA6hu0GNApDe+mDUAzmrrLZqaKBUS/nZLLP4xNQ5KhYTVx/Ix47MDKK3RiI5GREREZDYsnIg6wenCatz7xQFU1uuQFOGFJbMHwtVRJTqWWd0zqAuWzR4IDycVjpy/hCkf7UFmUbXoWERERERmwcKJqIOdLa7GPZ/vR3mtFglhnlg2ZyDcbKxoajKsux9WPTYUXXxdcOFSPaZ9vBfbTxeLjkVERER03Vg4EXWgrJIa3P3ZAZTWaBEb7IGv5gyCh5OD6FgdqluAG9Y8NhQDo3xQrdFjztJD2JjGiXtERERk3Vg4EXWQ82W1mPHZAZRUaxAT5I5vHhgETxfbLpqaeLuq8fXcQbitbyiMMvD8jydQWNkgOhYRERFRu7FwIuoAeeV1mPHZARRWNaB7gBu+fmAQvF3VomN1KrVKgX9OS0BimCeqGvR47ofjMBo5qpyIiIisk22eaEEkUH5FPe7+bD/yK+oR7e+Kbx4cBD83R9GxhHBQKvCf6X0w8YNd2H22FEv35mBO4+hyIiKyTA06A6oadKiq1zf+V4eqBj2q6nVwd1JheHd/uzsYSASwcCIyq7IaDWZ8th8XLtUj0tcF3z04GAHuTqJjCdXV3w0vTozF39ak4a2NGRjW3Q89At1FxyIishuyLOPEhUqkXay8rBiq/F1B9PtCSas3XvPxFBLQr4s3RvcKxJhegejq7wpJkjrppyESR5Jl2a56Z6qqquDp6YnKykp4eHiIjkM25sXVqfjmQC7CvJ3x/cPJCPFyFh3JIsiyjNlLD2H76RLEBntgzbyhUKvYKUxE1JFyy+qw+lg+1qTkI7u0tk33lSTAw8kBHs4q03+dHODupEJueR0yCi/faiLS1wWjewVidK8ADIj0gYOSr+9kPdpSG7BwIjKTs8XVGPfeLhiMMpY/NBiDo31FR7IoxVUNGPfeTlyq0+GREV3xwoQY0ZGIiGzOpVot1qUWYM2xfBw5f6n5eicHBZKjfeHr5nh5QeTsAA8nVeN/Tdd7OjvAVa2CQnHlVaQLl+qwNaMYm9OLsf9cGbSG/61QeTipMKJnAMb0CsDIHgF2MxSJrBcLp2tg4UQd5YFlh7A5vRhjegXi81n9RcexSBvTCvHI10cgScCKh5IxMMpHdCQiIqvXoDNgS3oxVh/Lx/bTxdA3DuJRSMDQbn6YmhSKm3oHdcgegjUaPXafKcHm9GJszShGea22+XtKhYT+XbwxNjYQo3sFIsrP1ezPT3S9WDhdAwsn6gj7s8pw16f7oVRI+PXp4egW4CY6ksX60w/H8cORCwj1csbGp2+Au43va0VE1BGMRhn7s8uw5lg+fkktRLVG3/y9uFAPTOkTilsTQxDg0Xnn2RqMMlLyLmFzejG2pBchs6jmsu9H+7tiTK9A3JoYgrhQz07LRXQtLJyugYUTmZvRKGPKx3tw4kIl7hkUgX9MjRcdyaLVaPSY8P5O5JXXY1rfMLxzZ6LoSEREVuN0YTVWH8vHTyn5KPjd/nihXs6YkhSCKX1C0d1CBvDkltVhS0YRtqQX40B2GXSG/62E/WlcDB4ZEc2hEiQcC6drYOFE5vZTSj6eWp4CV7US2/80Cv7u9jl6vC0O55Tjzk/2wSgDC+/piwnxwaIjERFZrLIaDVYevYDVxy4ivaCq+XoPJxUmJoRgalIo+nfxvuo5SZagqkGHXZmlWH0sH5vTiwAA43oH4t93JMKDnQckEAuna2DhROak0Rsw+p0duHCpHs+O7YEnRncXHclq/PvXDHy07Ry8XBzw69PDEdiJ7SRERNZAlmWsOJSHNzako6rB1IqnViowKsYfU5NCMSomAI4qpeCUbSPLMr49mItXfz4FrcGIKD9XLLq3H3oGWcYqGdkfFk7XwMKJzOmznVn4x4Z0BHo4Yvtzo+Cstq43MJG0eiNuW7gHaflVGN7DH8tmD2DLBhFRo3MlNfjLqlQcyC4HAPQK9sB9g7vg5vggeLlY/+azx/Mq8OjXR3CxsgHODkq8NS0ek/uEio5FdqgttQEH7RO1U0WdFv/degYA8OzYniya2kitUuC96X3gqFJgZ2YJvtp/XnQkIiLhtHoj/rvlDCa8vwsHssvh7KDEXyf2wtrHh2LGoAibKJoAIDHcC+uevAE3dPdDvc6Ap5an4JWfT7a4+S6RSCyciNrpv1vPoqpBj5ggd0zrFyY6jlXqFuCOBY37Of1jfTrOFte0cA8iItt1NPcSJv13N97ZlAmt3ogRPfzx2zPD8cAN0VDZ4KayPq5qLJ09EI+P6gYAWLo3B3d/th+Fvxt6QWRJbO9fIVEnyC2rw5f7cgAAL0yIgdKCT8i1dDOTI3FDdz9o9EY8syKFRxuJyO7UaPR4+ac0TFu4F6eLquHrqsb7d/XB0tkDEO7jIjpeh1IqJDw3ric+n9kf7k4qHDl/Cbf8dxf2nSsTHY3oD1g4EbXDv37NgM4gY1g3P4zo4S86jlVTKCS8fUcivFwckJpfiQ+2nBEdiYio02w+VYSx/9mBZfvOQ5aB2/uFYfP8EZjcJ9SuzvscExuItY8PQ0yQO0prtLj3iwP4dOc52Nmp+GThWDgRtVFKXgXWnSiAJAELbo6xqze2jhLo4YQ3Gve/+nj7WRw5Xy44ERFRxyquasC8b47igS8Po6CyAV18XfDNA4Pw9h2J8Ha1jfOY2irSzxWrHxuK25JCYTDKeGNDBh775ihqfre5L5FILJyI2kCWZbyxPh0AcFtSGHqHcOdzc7k5Phi3JYXCKAPPrDjON0oisklGo4zvDuZi9H92YH1qAZQKCY+M6IqNTw3H0G5+ouMJ56xW4p07E/H3KXFwUEr4Ja0Qt364G2eKqkVHI2LhRNQWv50qwsGccjiqFHhuXA/RcWzOK5N7I9TLGbnldfj72lOi4xARmdW5khrc9dl+LFiViuoGPRLCPPHz40PxwoQYTmb9HUmScN/gLljxcDKCPZ2QVVKLyR/twdrjF0VHIzvHwomolXQGI/75SwYAYO6wKAR7OgtOZHs8nBzwnzsTIUnAisN5+PVkoehIRETXTas34oMtZzDhvV04+LsR46seHcLOhWvoG+GNdU8Mw5CuvqjTGvDEd8fw2tpT0Bk4RIjEYOFE1ErLD+Yiq7QWPq5qPDKyq+g4NmtQtC8eGh4NAFiwKhXF1RxLS0TW68SFCtzy3134z6ZMaA22P2Lc3HzdHPHlnIF4tPF9d/GebMz4bD+Kq/jeQJ2P/2KJWqG6QYf3NpumvT09pjs8nBwEJ7Jt88f2QK9gD5TXavHnH09wqhIRWaW1xy/i9kX7kFlUY1cjxs1NpVTgz+Nj8Ml9/eDuqMKhnEuY+vFeVNRpRUcjO8PCiagVPtmRhbJaLaL9XHH3wAjRcWyeo0qJ96b3gVqlwLbTJfjmQK7oSERErSbLMt7ffAZPfHcMWr0Ro2MC7HLEuLmN6x2En58YhggfF+RX1GPBqlQeWKNOxcKJqAUFlfX4fHcWAOD58TFwYGtFp+gZ5I7nx/UEAPxjfTqyS2sFJyIialmDzoCnlqfg3c2ZAIAHb4jCpzP72+2IcXOL8nPFRzP6Nk/cW3EoT3QksiP8BEjUgnd+y0SDzoj+Xbwxrneg6Dh2Zc7QKAzp6ot6nQH/2ZQpOg4R0TWVVGsw47P9+Pn4RagUEt68LR4vToyFUsFVJnOKD/PEczeZDqy9uvYUzhbXCE5E9kJo4fTmm29iwIABcHd3R0BAAKZMmYLTp0+3eL8ffvgBMTExcHJyQnx8PDZs2NAJackenbpYhZVHLwAA/jKxF1ssOplCIeGvE2MBAOtOXOSbIxFZrNOF1Zjy0R4cza2Ah5MKX84ZyNbuDvTgDdEY2s10YO2p5ceg0RtERyI7ILRw2rFjB+bNm4f9+/dj06ZN0Ol0uOmmm1Bbe/WWnL179+Luu+/G3LlzcezYMUyZMgVTpkxBWlpaJyYne/HmL+mQZWBiQjD6RniLjmOXYkM8MDY2ELIMfLTtrOg4RER/sC2jGNMW7kV+RT2i/FyxZt5QDOFmth1KoZDwnzv7wNvFAScvVuHtX1s+8E50vSTZgs6qKykpQUBAAHbs2IHhw4df8TbTp09HbW0t1q1b13zd4MGD0adPHyxatKjF56iqqoKnpycqKyvh4eFhtuxke3ZmlmDm4oNwUErYMn8kInw5BUmU1AuVmPThbigkYOuzIxHp5yo6EhERZFnG0r05+Pu6UzDKwOBoHyy6tx+8XHg+U2fZfKoID3x5GACwbM5AjOjhLzgRWZu21AYWdY5TZWUlAMDHx+eqt9m3bx/GjBlz2XXjxo3Dvn37OjQb2ReDUcYbG9IBAPcNjmTRJFh8mCdG9fSHkatORGQhdAYj/vZTGl5dayqa7uwfhi/nDGLR1MnGxAZiZnIXAMCz3x9HaY1GcCKyZRZTOBmNRjz99NMYOnQo4uLirnq7wsJCBAZefoJ+YGAgCgsLr3h7jUaDqqqqyy5ELVl19AIyCqvh7qTCEzd2Ex2HADwxujsAYNWxfOSV1wlOQ0T2rLJehzlLD+Hr/bmQJOAvN8fgn9MSoFZZzMcqu/KXm3uhR6AbSms0+NMPxzminDqMxfwLnzdvHtLS0rB8+XKzPu6bb74JT0/P5kt4eLhZH59sT73WgLd/M/VKPz6qG0fIWoi+Ed64obsfDEYZH28/JzoOEdmp82W1uO3jPdh1phTODkp8cm8/PDS8K4cHCeTkoMQHdyc17/23bG+O6EhkoyyicHr88cexbt06bNu2DWFhYde8bVBQEIqKii67rqioCEFBQVe8/YIFC1BZWdl8ycvjvH+6ti92Z6GoSoNQL2fMGhIpOg79zpONq04/HslDfkW94DREZG8OZpdjykd7cK6kFsGeTvjhkWTc1PvKnz+oc8UEeeDFm3sBAN74JQPpBewwIvMTWjjJsozHH38cq1evxtatWxEVFdXifZKTk7Fly5bLrtu0aROSk5OveHtHR0d4eHhcdiG6mtIaDRbtaNrstiecHJSCE9HvDYj0QXK0L3QGGYu46kREnejHIxdwz+f7calOh4QwT/w0byjiQj1Fx6LfmZncBaNjAqDVG/Hkd8dQr+WIcjIvoYXTvHnz8PXXX+Pbb7+Fu7s7CgsLUVhYiPr6/x1JnjlzJhYsWND89VNPPYWNGzfinXfeQUZGBl555RUcPnwYjz/+uIgfgWzM+5vPoEajR3yoJyYlhIiOQ1fQtOq04lAeCisbBKchIltnNMr418YMPPfDcegMMm6OD8KKh5IR4OEkOhr9P5Ik4V+3J8Df3RFnimvwjw2nREciGyO0cFq4cCEqKysxcuRIBAcHN19WrFjRfJvc3FwUFBQ0fz1kyBB8++23+PTTT5GYmIgff/wRa9asueZACaLWOF9Wi28P5gIwnWiq4E7vFmlwtA8GRHpDazDik51cdSKijlOn1eOxb442n1f5+Khu+PDuvnBWsxvBUvm6OeI/dyYCAL7en4vfTl55eBhRe1jUPk6dgfs40dX8+9cMfLTtHG7o7oev5g4SHYeuYdeZEtz3xUE4qhTY9edRCHDnkV8iMq+Sag0eWHYIxy9UQq1U4K1p8bit77XPwybL8caGdHy6MwteLg7Y+NRwBHnyfYKuzGr3cSISxWiUsfpoPgBg+gBOXrR0w7r5ISnCCxq9EZ/vyhYdh4hsTFZJDaYt3IvjFyrh7eKAbx4cxKLJyjx3U0/EhXqgok6H+d+nwGC0q3UC6iAsnIgA7Msqw8XKBng4qTCmV2DLdyChJEnCkzeaznX6at95lHHDQyIykyPnyzFt4V7kltchwscFKx8dggGRPqJjURupVQp8cFcSnB2U2HuuDJ/uzBIdiWwACycimKYlAcCkxBBO0rMSI3v6Iz7UE/U6Az7fzVUnIrp+G9MKMeOzA7hUp0NimCdWPTYE0f5uomNRO0X7u+HVW3sDAN757TSO51WIDURWj4UT2b0ajR4b00wnj07rx1YMayFJUvOEvS/35qCiTis4ERFZs6V7svHoN0eg0RsxOiYA3z00GH5ujqJj0XW6o38YJsYHQ2+U8dTyY6jR6EVHIivGwons3obUAtTrDIj2d0VSuJfoONQGY3oFoFewB2q1BizmqhMRtYPRKOONDel4Ze0pyDJwz6AIfHJfP7ioVaKjkRlIkoQ3psYjxNMJOWV1eOXnk6IjkRVj4UR2r6lNb1rfMEgSR5BbE9O5Tt0AAEv25KCyXic4ERFZkwadAU8uP9Z8/svz43vi9SlxUCn58ciWeLo44L27kqCQTO/5Px+/KDoSWSm+MpBdyyuvw8HsckgScFvfUNFxqB3G9Q5Cj0A3VGv0WLY3R3QcIrISlXU6zFx8EOtOFMBBKeHd6Yl4bGQ3HkCzUQOjfPD4KNOBthdXpSKvvE5wIrJGLJzIrq08alptGtbND8GezoLTUHsoFBIeb5yw98XubFQ3cNWJiK7twqU6TFu0Fwezy+HuqMLS2QMxNYnnuNq6J0d3R98IL1Rr9Hh6RQr0BqPoSGRlWDiR3TIa5ebCaRr357BqE+ODEe3visp6Hb7cd150HCKyYGn5lZj68V6cLa5BkIcTfng0GUO7+YmORZ1ApVTg/buS4O6owpHzl/DfrWdFRyIrw8KJ7NahnHLkldfDzVGFcb2DRMeh66BUSHii8VynL3Zno5ZTk4joCnZklmD6J/tQUq1Bz0B3rJ43BDFBHqJjUScK93HB61PjAAALt59DUVWD4ERkTVg4kd1qWm2aGB8MZzX3brJ2kxJC0MXXBeW1WnxzgKtORHS57w/nYc7SQ6jVGjCkqy9+eDSZLdp2anKfUAyI9IbWYMTnu7gxLrUeCyeyS3VaPTakcu8mW6JSKjCv8cTfT3dmoV5rEJyIiCyBLMt4b3Mmnv/xBAxGGVP6hGDp7IHwcHIQHY0Eemyk6f3i2wO53AeQWo2FE9mlX08WokajR4SPCwZEeouOQ2YyNSkUYd7OKK3R4tuDuaLjEJFgOoMRf155Au9tPgMAeGxkV7w7vQ/UKn78sXcje/ojJsgdtVoDz42lVuMrB9mllUfyAXDvJlvjoFQ0H0X8ZMc5NOi46kRkr2o0ejyw7DC+P3wBCgl4fUocnh8fw9d8AmDaB/DRkV0BAEv2ZKNOy3NjqWUsnMjuXKyox55zpQC4d5MtmtYvFCGeTiiu1uD7w3mi4xCRAHnldbh94V7syCyBs4MSn97XH/cO7iI6FlmYifHBiPBxwaU6HVYc4vsFtYyFE9md1cfyIcvAoCgfhPu4iI5DZuaoUjYfRVy4/Rw0eq46EdmT/VllmPzRHmQUVsPPzRHfPTQYY2IDRcciC6RSKvDwiGgAwGc7s6DVc18nujYWTmRXZFnGyiOmaXq3cyiEzbqjfzgCPRxRUNnQ3JZJRLbvmwPnce/nB1Beq0V8qCfWPjEUfcK9RMciCzatbxj83R1xsbIBP6Xw/YKujYUT2ZVjeRXIKq2Fs4MSE+KDRcehDuLkoMTDw02rTh9tOwsdd4cnsmk6gxF/W5OGF1enQW+UMSkxBN8/zHHj1DInByXmDosCACzacQ5Goyw4EVkyFk5kV35sXG2aEBcEN0eV4DTUke4eGAE/N0fkV9Rj9VEeRSSyVeW1Wtz3xQF8tf88JAn407ie+OCuPtyfj1rtnkERcHdS4VxJLX47VSQ6DlkwFk5kNxp0Bqw7fhEA2/TsgbNaiYeHm3rXP9p+FnquOhHZnNOF1Zj80W7szyqHq9o0BGLeqG6cnEdt4u7kgJnJpuEhC7efhSxz1YmujIUT2Y3N6UWoatAj1MsZg6N9RcehTnDP4Aj4uKpxvqwOPzcWzURkG347WYjbPt6DvPJ6RPi4YPW8oRjLIRDUTrOHRsFRpcDxC5XYd65MdByyUCycyG40telNTQqFQsGjkfbARa3CAzeYetc/3HoWBvauE1k9WZbx4dYzeOirI6jVGpAc7Yuf5g1Fj0B30dHIivm5OWL6gHAAwMfbzwlOQ5aKhRPZheKqBuzMLAEATGObnl2ZmRwJLxcHZJXWYmNaoeg4RHQd6rUGPPHdMbz9WyYAYFZyF3w5dyC8XdWCk5EtePCGaCgVEnafLcWJCxWi45AFYuFEdmFNSj6MMtCvizei/FxFx6FO5OaowszGjS+X7c0RG4aI2u1iRT1uX7QX604UQKWQ8MbUeLw6OQ4OSn6UIfMI93HB5MQQAKZ9AIn+P77akM2TZbm5TW9aX6422aMZg7pAqZBwMKcc6QVVouMQURsdOV+OWz/cjZMXq+DjqsY3DwzCjEERomORDXqkcQP1jScLcba4RnAasjQsnMjmpeVXIbOoBo4qBSYmcO8mexTk6YTxvYMAAF/uOy84DRG1xfeH8nDXp/tRWqNFTJA7fpo3FIM44Ic6SI9Ad4zpFQhZBj7ZwVUnuhwLJ7J5K4+aVptu6h0ET2cHwWlIlKZRs2uO5aOyXic4DRG1RG8w4rW1p/D8yhPQGWSM7x2ElY8OQbiPi+hoZOMeG2VadVqTko+LFfWC05AlYeFENk2rN+KnFNPmp9P6hgpOQyINjPJBz0B31OsMza2bRGSZKut0mL30EBbvyQYAPDW6Oz6+py9cuXE5dYK+Ed4YFOUDnUHG57uyRcchC8LCiWza1oxiXKrTIdDDETd09xcdhwSSJAkzh5hWnb7alwMjR5MTWaTzZbWYunAPdp0phbODEgvv6YtnxvbgNhLUqR4b1Q0A8N3BXJTXagWnIUvBwolsWlOb3pSkUCj5pmv3pvQJhbuTCjllddh5pkR0HCL6fw7nlGPqx3uRVVKLYE8n/PhoMibE89xU6nzDu/uhd4gH6nUGTmSlZiycyGaV1WiwLaMYAHA7p+kRAFdHFW5v3MfrKw6JILIoP6XkY8bnB1Beq0VcqAfWzBuK3iGeomORnZIkCY82TthbujcHtRq94ERkCVg4kc36KeUi9EYZiWGe6M4d5anRfY17Om09XYzcsjrBaYhIlmX8d8sZPLU8BVq9EWNjA/H9w8kI9HASHY3s3IS4YET5uaKyXofvDuaKjkMWgIUT2aymNr1p/bjaRP8T7e+G4T38IcvA1we46kQkklZvxHM/nMA7mzIBAA8Mi8Kie/vBRc0hECSeUiHh4eHRAIDPdmVBozcITkSisXAim5ReUIWTF6vgoJQwKSFEdByyMDMbV51WHMpDvZZvhEQiVNRpcd8XB7Dy6AUoFRJenxKHv94Sy/NRyaJM7RuKQA9HFFVpsOZYvug4JBgLJ7JJKxvHTY+OCYS3q1pwGrI0o2ICEObtjMp6HdYevyg6DpHdySmtxW0f78WB7HK4Oaqw+P4BuLfxgAaRJXFUKfHAMNOq06IdWTBwIqtdY+FENkdnMGJNiunD8O1s06MrUCqk5nOdlu3LgSzzjZCosxzKKcfUj/cgq7QWoV7O+PHRZIzowe0iyHLdPSgCns4OyC6txa8nC0XHIYFYOJHN2ZlZgtIaDXxd1RjRk2/GdGV39g+Ho0qBkxercDS3QnQcIrvwU0o+7vnsAC7V6ZAQ5onV84YgJshDdCyia3JzVGFWsulg28fbz/Jgmx1j4UQ2p2koxOQ+oXBQ8lecrszbVY1bE03nv325L0dsGCIbJ8sy3t/cODnPYMS43oFY8VAyAtw5OY+sw/1Do+DsoERafhV2nSkVHYcE4adKsikVdVpsPtW4dxPb9KgFs4ZEAgA2pBagpFojNgyRjdLoDXj2++N4d7Npct5Dw6Ox8J5+cFYrBScjaj0fVzXuGhgOAFi4/ZzgNCQKCyeyKWuPX4TWYESvYA/EhrD9g64tLtQTSRFe0BlkLOceHURmd6lWi/u+OIhVx/KhVEh4Y2o8/nJzLyg4OY+s0AM3REOlkLAvqwzHci+JjkMCsHAim/LjUdOo0Gl9QwUnIWsxKzkSAPDNgVzoDEaxYYhsSHZpLW5buBcHs8vh7qjCkvsHYMagCNGxiNot1MsZU5JMny8+5qqTXWLhRDbjbHENjudVQKWQml/YiFoyIT4Ifm5qFFY1YNOpItFxiGzCwWzT5Lzs5sl5QzCck/PIBjwyIhqSBGw6VYQzRdWi41AnY+FENqNpKMTInv7wc3MUnIashaNKibsHmo6Cc0gE0fXbdKoI935+ABV1OiQ2Ts7rGeQuOhaRWXQLcMdNsYEAgIU7uOpkb1g4kU0wGGWsaiycpvXlUAhqmxmDIqBUSNifVY7ThTyCSNReh3LK8fi3R6E1GDG+dxCWc3Ie2aBHR3YDAPycchEXLtUJTkOdiYUT2YSD2eUoqtLA09kBN/YKEB2HrEywp3PzEUSuOhG1z5miasxdeggavRFjegXiwxlJnJxHNqlPuBeGdPWF3ijjs51ZouNQJ2LhRDah6dyUsbGBcFTxjZrabmbjkIjVx/JR1aATG4bIyhRU1mPW4oOoatCjb4QX/nt3ElTcR49s2KMjuwIAVh7NR4POIDgNdRa+qpHVk2UZm9ILAZgKJ6L2GBztgx6BbqjTGrDyyAXRcYisRmW9DvcvPoSLlQ3o6u+KL2YN4EoT2byhXf0Q6uWMGo0eWzOKRcehTsLCiaxeZlEN8srr4ahS4IbufqLjkJWSJAn3Na46fbXvPIxGWWwgIivQoDPgwS8P43RRNQLcHbFszkB4u6pFxyLqcAqFhFsSgwGYznUi+8DCiazeplOm1aZh3fzgolYJTkPW7LakULg7qpBVWos950pFxyGyaAajjGdWpDTv07RszkCEebuIjkXUaSYnmrY+2Xq6mC3edoKFE1m9TemmJfIxbNOj6+TqqMK0fqapjMv2nhechshyybKMV9eexC9phVArFfhkZj/0CvYQHYuoU/UKdke3ADdo9Ub8mlYoOg51AhZOZNWKqhpwPK8CADCa0/TIDO4d3AUAsCWjCHnlHDNLdCUfbz+HL/edhyQB/5meiCFd2SZN9keSJExODAEA/Hyc7Xr2gIUTWbUtjatNfcK9uFcImUW3ADfc0N0Psgx8cyBXdBwii/PD4Tz8+9fTAIC/TYzFLQkhghMRiTOpsXDac7YUJdUawWmoo7FwIqvWdH4Tp+mROd3XuOq04lAux8wS/c6208V4YVUqAODhEdGYMyxKcCIisSL9XJEY5gmjDGxILRAdhzoYCyeyWrUaPfacKwPAwonMa3SvQIR6OeNSnQ5r2X5BBABIyavAY18fhcEoY2pSKP48LkZ0JCKLcGsf05CIn1LyBSehjsbCiazWrjMl0OqN6OLrgu4BbqLjkA1RKqTmc52+3HcesszR5GTfsktrMWfpIdTrDLihux/+OS0BCoUkOhaRRbglIRiSBBzNreC5sTaOhRNZrU2nGqfp9QqEJPENnMxr+oBwqFUKpOZXIqVxAAmRPSqubsDMxQdQXqtFfKgnFt7bD2oVPz4QNQn0cEJytC8ADomwdXzlI6ukNxixNaMIANv0qGP4uKoxqfGk9y/3cTQ52acajR5zlh5CXnk9InxcsPj+AXBz5H55RP/frY1DItjebdtYOJFVOppbgUt1Oni5OKB/F2/RcchGzRpiatdbf6IApTWclkT2Ras34tGvjyAtvwq+rmp8OWcg/N0dRcciskgT4oLhoJSQUViN04XVouNQB2HhRFapaZrejT0DoFLy15g6RkKYFxLDvaA1GLHiUJ7oOESdxmiU8fyPx7HrTClc1Eosvn8AIv1cRccislieLg4Y0cO0n+TPxzkkwlbxEydZHVmWsemUqU1vDNv0qIPNSjatOn29/zz0BqPgNESd462NGViTchEqhYSP7+mLxHAv0ZGILN6tff63GS6HCtkmFk5kdc6V1CCnrA5qpQLDe/iLjkM27ub4YPi6qlFQ2YDN6UWi4xB1uM93ZeHTnVkAgH9OS8DIngGCExFZhzG9AuCiViKvvB7HOFTIJrFwIqvzW+Nq05BuvjxJmTqck4MSdw0MB8AhEWT7tmYU4fX16QCA58f3xLR+YYITEVkPF7WqeWDVzykcEmGLWDiR1dnc1KbXi2161DlmDOoChQTsPVeGzCKe9Eu2qVajx4ur0wAA9wyKwKMjugpORGR9Jje26607UcD2bhvEwomsSkm1pnn5m4UTdZZQL2fcFBsEAPhkR5bgNEQd4/0tZ1BQ2YAwb2f8dWIs98cjaodh3fzh5eKA0hoN9meVi45DZsbCiazK1owiyDKQEOaJIE8n0XHIjjwy0nT0/aeUfORX1AtOQ2Re6QVV+GJ3NgDg75Pj4KxWCk5EZJ3UKgVujg8GwOl6tkho4bRz505MmjQJISEhkCQJa9asuebtt2/fDkmS/nApLCzsnMAkXNM0vbFcbaJO1ifcC0O6+kJvlPHZTq46ke0wGmW8uDoVBqOMCXFBGBXDYRBE16NpM9xf0gqh0RsEpyFzElo41dbWIjExER999FGb7nf69GkUFBQ0XwIC+CJvD+q1Buw6UwqAY8hJjMdGdgMALD+UizJuiEs2YsXhPBzNrYCrWomXJsWKjkNk9QZG+iDIwwnVDXpsP10iOg6ZkdDCacKECXj99dcxderUNt0vICAAQUFBzReFgh2H9mDXmRJo9EaEeTsjJshddByyQ0O7+SIhzBMNOiOW7MkRHYfoupXWaPDWLxkAgGfG9kCwp7PgRETWT6GQMCmxsV2P0/VsilVWHH369EFwcDDGjh2LPXv2XPO2Go0GVVVVl13IOjXtoTOmVyBPWiYhJEnCY43nOi3bl4PqBp3gRETX540N6ais16FXsAfuHxIpOg6Rzbg1MRSA6bNLjUYvOA2Zi1UVTsHBwVi0aBFWrlyJlStXIjw8HCNHjsTRo0evep8333wTnp6ezZfw8PBOTEzmYjDK2JJeDAC4iW16JNBNsUHo6u+K6gY9vj2QKzoOUbvtO1eGVUfzIUnAG1PjoFJa1UcCIosWF+qBaD9XaPRG/HaS5+LbCqt6lezZsycefvhh9OvXD0OGDMHixYsxZMgQvPvuu1e9z4IFC1BZWdl8ycvL68TEZC4peZdQVquFh5MKA6J8RMchO6ZQSHikcX+bz3dno0HHE3/J+mj1Rvx1TSoAYMbACCRFeAtORGRbJEnCpMYhET8fZ7uerbCqwulKBg4ciLNnz171+46OjvDw8LjsQtbnt8ZpeqNiAuDAo6Ik2OQ+oQjxdEJJtQYrj14QHYeozT7blYVzJbXwc1Pj+XExouMQ2aRbGzfD3XWmlAOFbITVfwJNSUlBcHCw6BjUwTaf+t/5TUSiqVUKPDg8GoBpQ1zuDk/WJLesDh9sOQMA+OvEWHi6OAhORGSbuvq7IS7UAwajjA1pbNezBUILp5qaGqSkpCAlJQUAkJ2djZSUFOTmms4bWLBgAWbOnNl8+/feew8//fQTzp49i7S0NDz99NPYunUr5s2bJyI+dZKskhqcK6mFg1LCiJ7+ouMQAQDuGhABH1c1csvrsD61QHQcolaRZRl/+ykNGr0RQ7r6YnLjEXEi6hhNezqt5XQ9myC0cDp8+DCSkpKQlJQEAJg/fz6SkpLw0ksvAQAKCgqaiygA0Gq1ePbZZxEfH48RI0bg+PHj2Lx5M0aPHi0kP3WOpml6g6N94eHEI6NkGZzVyuYpZAu3n4Msy2IDEbXCL2mF2JFZArVSgb9PieOEUqIONikxBJIEHMwpR35Fveg4dJ1UIp985MiR1/ywsXTp0su+fv755/H88893cCqyNJsa2/TGcpoeWZhZyZH4ZMc5ZBRWY9vpYtwYw99Rslw1Gj1eXXsSAPDIiGh09XcTnIjI9gV7OmNApA8OZpdj3fGLeLhxuBBZJ6s/x4lsW1mNBkfOXwIAjOb5TWRhPF0ccM/gLgCAj7edE5yG6Nr+81smiqo06OLrgsdGdRMdh8huNLXEcrqe9WPhRBZta0YxjDLQO8QDoV7c0Z4sz9xhUVArFTh8/hIOZpeLjkN0RWn5lVi6NxsA8NrkODg5KAUnIrIfN8cFQ6WQcPJiFc4W14iOQ9eBhRNZtE2cpkcWLtDDCdP6hQEAPt5+9a0RiEQxGGW8uDoVRhm4JSEYI3pwyA5RZ/J2VWN44787rjpZNxZOZLEadAbsOlMKgOc3kWV7ZEQ0FBKw/XQJTl6sFB2H6DLfHszF8QuVcHdU4W+3xIqOQ2SXmqbr/ZySz2FCVoyFE1msPWdLUa8zIMTTCb1DuHExWa4uvq6YmGB6U1y4nec6keUorm7AvzZmAACeG9cTgR5OghMR2aexsYFwclAgp6wOqfk8wGatWDiRxWoaQz4mNpAjc8niPdo4KWlDagFySmsFpyEyeWN9Oqob9IgP9cS9jYNMiKjzuTqqmk87+Il7OlktFk5kkYxGGZvTiwGwTY+sQ2yIB0b19IdRBj7ZmSU6DhH2nC3FmpSLUEjAG1PjoVTwABSRSE3teutOXITByHY9a8TCiSzS8QsVKKnWwN1RhUFRvqLjELXKoyNNI55XHrmAoqoGwWnInjXoDPjrmjQAwH2DuyA+zFNwIiIa0dMfHk4qFFVpOIXVSrFwIovUNE1vRE9/qFX8NSXrMDDKB/27eENrMOKL3dmi45Ad+2RHFrJLa+Hv7ohnx/UUHYeIADiqlJgQFwwA+Pl4vuA01B78REoWqen8JrbpkbV5bJTpXKdv9p9HZZ1OcBqyR9mltfiocTT+S7fEwsPJQXAiImrStBnuhtRCaPVGwWmorVg4kcU5X1aLzKIaqBQSRvYIEB2HqE1G9QxATJA7arUGLNuXIzoO2RlZlvHST2nQ6o24obsfbkkIFh2JiH5nULQvAtwdUVmvw87MEtFxqI1YOJHFaWrTGxjlA08XHikl6yJJEh4daVp1WrInG3VaveBEZE/WnSjArjOlUKsU+PvkOE4kJbIwSoWEWxq3r+BmuNaHhRNZnKbCiW16ZK0mxgcjwscFl+p0WH4wT3QcshNVDTq8tu4UAODxUd0Q6ecqOBERXcmtje16m04V8eCalWHhRBblUq0Wh89fAoDm/Q6IrI1KqcDDI6IBAJ/tymIfO3WKpXtyUFKtQbSfa/PvHxFZnsQwT3TxdUG9ztB8sJisAwsnsijbThfDYJQRE+SOcB8X0XGI2m1a3zD4uzuioLIBa1I4PYk6llZvxFf7zwMAnhrTHY4qpeBERHQ1kiQ17+m0lu16VoWFE1kUTtMjW+HkoMQDw6IAAIt2nONmh9ShNqQWoKRagwB3x+Zxx0RkuZoKpx2ZJaio0wpOQ63FwokshkZvwI7TpgkzLJzIFswYFAEPJxWySmqx6VSh6Dhko2RZxpI9pn3D7hvchXvfEVmB7oHu6BXsAZ1Bxi9pfH+wFtf16qrVanH69Gno9Tyxja7f3nNlqNUaEOjhiLgQ7nJP1s/dyQEzkyMBAB9vPwdZ5qoTmd/R3Aocv1AJtUqBGYMiRMcholZqWnX6OYXtetaiXYVTXV0d5s6dCxcXF/Tu3Ru5ubkAgCeeeAJvvfWWWQOS/djceILkmF6BUCg4Qpdsw+yhkXByUODEhUrsOVsmOg7ZoKbVpsmJIfB1cxSchohaa2K8qa32YE45N0y3Eu0qnBYsWIDjx49j+/btcHJyar5+zJgxWLFihdnCkf0wGuXm85vGsE2PbIivmyPuGmBaBfh4+1nBacjWFFTWN7f5zB4aJTgNEbVFhK8Lega6w2CUsT2zWHQcaoV2FU5r1qzBhx9+iGHDhl22uV7v3r1x7tw5s4Uj+5F2sRJFVRq4qpUY0tVXdBwis3pweDRUCgl7z5UhJa9CdByyIV/tOw+DUcagKB/EhniIjkNEbTS6VwAAcCy5lWhX4VRSUoKAgIA/XF9bW8tdyqldml4whvfw5xhdsjmhXs6Y3CcUAPDxNq46kXk06Az47qCpVZ6rTUTWqanLZsfpEu75ZwXaVTj1798f69evb/66qVj6/PPPkZycbJ5kZFeaCidO0yNb9ejIaEgS8NupIqTlV4qOQzZgzbF8XKrTIczbma+dRFaqT5gX/NzUqNbocSinXHQcaoGqPXd64403MGHCBJw6dQp6vR7vv/8+Tp06hb1792LHjh3mzkg2Lq+8DhmF1VAqJNwY88eVTCJb0C3AHbckhGDt8Yt4YdUJrHlsKFRKjo2m9jGNIM8BAMxKjoSSA3WIrJKi8bPP94cvYNOpIgzt5ic6El1Du961hw0bhpSUFOj1esTHx+O3335DQEAA9u3bh379+pk7I9m4DakFAICBkT7wclELTkPUcf52Sy94OKmQll+FL3Zni45DVmzfuTKcLqqGi1qJOweEi45DRNdhTC/TivHm9CJuW2Hh2rXiBABdu3bFZ599Zs4sZKc2NE6EujmBu92TbQtwd8Jfb4nF8z+ewH82ZWJc7yBE+rmKjkVWaHHjatO0vmHwdHYQG4aIrsuw7n5wVClw4VI9Motq0DPIXXQkuop2rTht2LABv/766x+u//XXX/HLL79cdyiyHxcu1eF4XgUkCRjfO0h0HKIOd0e/MAzr5geN3ogXVp3g0UVqs/NltdiSYTov9P6hkWLDENF1c1Grmlv0mrZmIcvUrsLphRdegMFg+MP1sizjhRdeuO5QZD82Nq42DYz0gb87N24k2ydJEt6YGg9nByX2Z5Vj+aE80ZHIyizbex6yDIzo4Y+u/m6i4xCRGTS163EsuWVrV+F05swZxMbG/uH6mJgYnD3LUbvUek3nN90czzY9sh8Rvi549qYeAIA3NqSjqKpBcCKyFjUaPX44bCq2Z3O1ichmNO3ndPxCBYqr+Z5gqdpVOHl6eiIrK+sP1589exauruzXp9a5WFGPo7mNbXpxbNMj+zJ7aBQSw71Q3aDH39aksWWPWuXHw3mo1ugR7e+K4d39RcchIjMJ9HBCQpgnZBnYllEsOg5dRbsKp8mTJ+Ppp5/GuXPnmq87e/Ysnn32Wdx6661mC0e2ralNr38XbwR6OAlOQ9S5lAoJ/5wWD5VCwm+nivBL478HoqsxGmUs23ceADB7SCQUHEFOZFP+167HwslStatw+te//gVXV1fExMQgKioKUVFR6NWrF3x9ffH222+bOyPZqF/S2KZH9i0myAOPjewKAHjppzRU1GkFJyJLtj2zGNmltXB3UuG2vmGi4xCRmTUVTrvPlqBB98dZAiReu8aRe3p6Yu/evdi0aROOHz8OZ2dnJCQkYPjw4ebORzaqqKoBh89fAsA2PbJv827shg1phThbXIPX16fj7TsSRUciC9W04e1dA8Lh6tju3USIyEL1CnZHiKcTLlY2YM/ZUoxuLKTIcrR723pJknDTTTfhT3/6Ex5//HEWTdQmG9MKIctA3wgvBHs6i45DJIyjSol/TouHJAE/HrmAXWdKREciC3SmqBq7zpRCIQEzkyNFxyGiDiBJEsbE/m8zXLI87T5ktWXLFmzZsgXFxcUwGo2XfW/x4sXXHYxs23pO0yNq1q+LD2YlR2Lp3hwsWJWK354ZDhc1VxTof5bszQEAjI0NRLiPi9gwRNRhxvQKxJf7zmNLejGMRpnnMlqYdq04vfrqq7jpppuwZcsWlJaW4tKlS5ddiK6luLoBh3LKAQATWDgRAQD+NK4nQr2cceFSPd75LVN0HLIgFXVarDp6AYBpGiMR2a5B0T5wVStRXK1Ban6l6Dj0/7TrkOaiRYuwdOlS3HfffebOQ3bg15NFkGWgT7gXQr3YpkcEAK6OKvxjahzuX3IIS/Zk45aEYCRFeIuORRZg+aE8NOiM6BXsgUFRPqLjEFEHclQpMaKnPzakFmJzehESw71ER6LfadeKk1arxZAhQ8ydhezEhhNNbXocCkH0eyN7BmBqUiiMMvDCylRo9caW70Q2TW8w4svGNr3ZQyMhSWzbIbJ1TdP1NqdzLLmlaVfh9MADD+Dbb781dxayA6U1GhzILgMATIhjmx7R//e3W2Lh66rG6aJqLNx+ruU7kE377VQRLlY2wNdVjVsTQ0THIaJOMKpnABQSkF5QhQuX6kTHod9pV6teQ0MDPv30U2zevBkJCQlwcHC47Pv/+c9/zBKObM+vJwthlIGEME+e4Ex0BT6uarx8a288+d0xfLjtDG6OD0L3QHfRsUiQJXuyAQAzBkXAyUEpOA0RdQZvVzX6d/HBwZxybEkvxqwhkaIjUaN2rTidOHECffr0gUKhQFpaGo4dO9Z8SUlJMXNEsiW/pBYC4GoT0bVMSgjG6JgA6Awy/rzyBAxGWXQkEiAtvxKHci5BpZBw7+AuouMQUScaExsAgGPJLU27Vpy2bdtm7hxkB8prtdiXZWrT4/lNRFcnSRJenxqHA//ZiaO5FfhqXw7u5zQ1u7O4cbVpYkIwAj2cBKchos40ulcg3tiQgf1ZZahu0MHdyaHlO1GHa/cGuERt9dvJQhiMMnqHeKCLr6voOEQWLdjTGX+eEAMA+Nevp9nnbmdKqjVYd9w0SIcjyInsT1d/N0T7uUJnkLEzs1R0HGrU7h0WDx8+jO+//x65ubnQarWXfW/VqlXXHYxsz4Y0U5seN70lap17BkZgbcpFHMwpx19Wp2HZ7AGcqmYnvjlwHlqDEUkRXujDccREdmlMbCA+3ZmFLelFmJjAz06WoF0rTsuXL8eQIUOQnp6O1atXQ6fT4eTJk9i6dSs8PT3NnZFsQEWdFnvPmo6YTIhjmx5RaygUEt6aFg+1SoGdmSVYfSxfdCTqBBq9AV/vzwXA1SYiezY6xnSe09bTxdAbuD2FJWhX4fTGG2/g3Xffxdq1a6FWq/H+++8jIyMDd955JyIiIsydkWzAb6eKoDfKiAlyR7S/m+g4RFYj2t8NT43uDgB4bd0plNZoBCeijrb+RAFKazQI8nDigSYiO9avize8XBxQUafDkfOXRMchtLNwOnfuHCZOnAgAUKvVqK2thSRJeOaZZ/Dpp5+aNSDZhl9STb36E9mmR9RmDw2PRmywByrqdHh17SnRcagDybKMJXtyAAD3JXeBg5KnIhPZK5VSgRt7mladtmRwM1xL0K5XZG9vb1RXVwMAQkNDkZaWBgCoqKhAXR1PYKbLVdbrsLupTY+FE1GbOSgV+Oe0BCgkYO3xi9h8iuNpbdWR85eQml8JR5UCdw9kBweRvRvdKxAA+LpvIdpVOA0fPhybNm0CANxxxx146qmn8OCDD+Luu+/G6NGjzRqQrN/mU0XQGWT0CHRDtwC26RG1R3yYJx4cHg0A+OuaNJTXalu4B1mjptWmKX1C4eOqFhuGiIQb3sMPDkoJWaW1OFdSIzqO3WtX4fThhx/irrvuAgC8+OKLmD9/PoqKijBt2jR88cUXZg1I1u+XNFObHqfpEV2fZ8b0QJSfKwqrGvDIV0eg0RtERyIzulhRj40nTdNHZw+LFBuGiCyCu5MDBkf7AgC2cDNc4dpVOPn4+CAkJMT0AAoFXnjhBfz8889455134O3tbdaAZN2qGnTN+w+wcCK6Pk4OSnxyXz+4O6pwMKccC1amQpZl0bHITL7cdx4Go4zkaF/EBHmIjkNEFmJMc7sez3MSrV2Fk1KpRHHxH//nlZWVQalUXncosh1b04uhNRjRLcANPQLdRcchsno9At3x0T19oVRIWHUsHx9uPSs6EplBvdaA7w42jSCPFBuGiCzK6F6mARGHz5fjEtu0hWpX4XS1I5wajQZqNXuy6X/WN07Tu5kjdYnMZngPf7w2uTcA4J1Nmfj5+EXBieh6bUgtQGW9DuE+zs0ngxMRAUCYtwt6BXvAKAPbTnPVSSRVW278wQcfAAAkScLnn38ON7f/nehvMBiwc+dOxMTEmDchWa0ajR47MksAcJoekbndM6gLckpr8dmubDz3w3GEejmjXxe2SlurpuL3jn7hUCokwWmIyNKM6RWA9IIqbE4vwm19w0THsVttKpzeffddAKYVp0WLFl3WlqdWqxEZGYlFixaZNyFZra0ZxdDqjYj2c0VMENv0iMzthQm9kFNWh02nivDQl4exZt5QhPu4iI5FbVRWo2nesmFSYojgNERkicb0CsR/t57FzsxSaPQGOKp4aowIbSqcsrOzAQCjRo3CqlWrOAiCrmnDCVOb3oT4IEgSj6ASmZtSIeH9u/rgjkX7cPJiFWYvPYSVjw6Bp7OD6GjUBhvSCmEwyogP9USUn6voOERkgeJDPRHg7ojiag0OZJVjeA9/0ZHsUrvOcdq2bdtlRZPBYEBKSgouXbpktmBk3Wo1+uY+XE7TI+o4LmoVvpg1AEEeTjhbXIN53xyFzmAUHYvaYG1jm96tXG0ioqtQKKTmIRGbOZZcmHYVTk8//XTzfk0GgwHDhw9H3759ER4eju3bt5szH1mpbaeLodEb0cXXBbHBHKtL1JGCPJ3w+az+cFErsftsKV76KY1jyq1EQWU9DuWUAwAmJvAgExFdXdNY8i3pxXyNF6RdhdMPP/yAxMREAMDatWuRk5ODjIwMPPPMM3jxxRfNGpCs0y+ppk0cJ8QFs02PqBPEhXrig7uSIEnAdwfz8PmubNGRqBXWHS+ALAMDI30Q4uUsOg4RWbCh3fzg5KBAfkU90guqRcexS+0qnMrKyhAUZBovvWHDBtxxxx3o0aMH5syZg9TUVLMGJOtTrzVga4apTW8i2/SIOs2Y2ED8dWIsAOCNX9Lx68lCwYmoJWtPmNr0JiXytZKIrs3JQYlh3UznNrFdT4x2FU6BgYE4deoUDAYDNm7ciLFjxwIA6urquAEuYfvpYtTrDAjzdkZcKNv0iDrTnKGRuHdwBGQZeHp5ClIvVIqORFeRXVqLExcqoVRIPBeUiFplbKzpPKctLJyEaFfhNHv2bNx5552Ii4uDJEkYM2YMAODAgQPcx4mwIc10lHtiPNv0iDqbJEl4ZVJv3NDdD/U6A+YuO4SCynrRsegK1jUOhRjazQ++bo6C0xCRNbgxJhCSBBy/UImiqgbRcexOuwqnV155BZ9//jkeeugh7NmzB46Ophd8pVKJF154wawBybo06AzNR0G46S2RGCqlAh/d0xc9At1QXK3BnKWHUavRi45FvyPLcvOmt5M4FIKIWsnf3RGJYV4ATEMiqHO1q3ACgNtvvx3PPPMMwsL+t3vxrFmzMHnyZLMEI+u0I7MEdVoDQr2ckRjmKToOkd3ycHLAF7MGwM9NjfSCKjz53TEYjJzCZCkyCqtxprgGapUC4+KCRMchIisyNrZpuh7b9TpbqwunDz74AA0NDc1/vtaltXbu3IlJkyYhJCQEkiRhzZo1Ld5n+/bt6Nu3LxwdHdGtWzcsXbq01c9HHe+X1MZNb+O46S2RaOE+Lvh0Zn84qhTYklGMf6xPFx2JGjXt3TSqpz88nLhhMRG1XtNY8t1nS1GnZTdBZ1K19obvvvsu7rnnHjg5OeHdd9+96u0kScKTTz7Zqsesra1FYmIi5syZg9tuu63F22dnZ2PixIl45JFH8M0332DLli144IEHEBwcjHHjxrX2R6EO0qAzYHPjsjHb9IgsQ98Ib7xzZyIe//YYFu/JRpSfC+5LjhQdy67Jsvy7aXrc9JaI2qZHoBvCvJ1x4VI9dp8pxU29uWrdWVpdOGVnZ1/xz9djwoQJmDBhQqtvv2jRIkRFReGdd94BAPTq1Qu7d+/Gu+++y8LJAuw+U4oajR7Bnk5ICvcSHYeIGt2SEILzZXX496+n8craUwj3ccHIngGiY9mtY3kVyCuvh6taidExgaLjEJGVkSQJY3oFYuneHGxJL2bh1IlaXTjNnz+/VbeTJKm5sDG3ffv2NU/wazJu3Dg8/fTTV72PRqOBRqNp/rqqqqpDshGwIc3Upjc+LggKBdv0iCzJYyO7Iru0Fj8euYDHvz2GlY8OQc8gd9Gx7FJTm97Y2EA4q7mFBxG13djYxsIpowhGo8zPXZ2k1YXTsWPHLvv66NGj0Ov16NmzJwAgMzMTSqUS/fr1M2/C3yksLERg4OVH5wIDA1FVVYX6+no4O/9x1/U333wTr776aodlIhON3oBNp0wnKXI/EiLLI0kS3pgaj7zyOhzILsecpYewet4QBLg7iY5mVwxGGetOmA4ysU2PiNprQKQP3B1VKK3RIuVCBfpGeIuOZBdaPRxi27ZtzZdJkyZhxIgRuHDhAo4ePYqjR48iLy8Po0aNwsSJEzsyb5stWLAAlZWVzZe8vDzRkWzS3rNlqG7QI8DdEf34j5fIIqlVCnxyXz9E+bkiv6IeD391BA06g+hYduVAVhlKqjXwdHbADd39RcchIiulVikwoqfpNYTT9TpPu8aRv/POO3jzzTfh7f2/D8je3t54/fXXO6xNDwCCgoJQVHT5L0dRURE8PDyuuNoEAI6OjvDw8LjsQua3/nfT9LhcTGS5vFzU+GJWf3g6O+BYbgWe//EEZJljyjtL01CIm+ODoFa1e0cQIqLmseSbT3E/p87SrlftqqoqlJSU/OH6kpISVFdXX3eoq0lOTsaWLVsuu27Tpk1ITk7usOeklmn1Rvx2shAA2/SIrEG0vxsW3tsXKoWEn49fxAdbzoqOZBe0eiM2pJpeKyclsE2PiK7PyB4BUCoknC6qRl55neg4dqFdhdPUqVMxe/ZsrFq1ChcuXMCFCxewcuVKzJ07t1VjxZvU1NQgJSUFKSkpAEzT+lJSUpCbmwvA1GY3c+bM5ts/8sgjyMrKwvPPP4+MjAx8/PHH+P777/HMM8+058cgM9mXVYaqBj383BzRP9JHdBwiaoUhXf3w+pQ4AMC7mzObBxZQx9l9tgSV9Tr4uztiULSv6DhEZOU8XRwwINLU/bWZ7Xqdol2F06JFizBhwgTMmDEDXbp0QZcuXTBjxgyMHz8eH3/8casf5/Dhw0hKSkJSUhIA0+S+pKQkvPTSSwCAgoKC5iIKAKKiorB+/Xps2rQJiYmJeOedd/D5559zFLlgG040TdMLhJJtekRW466BEXhgWBQA4LkfjuNY7iXBiWzbzymm4nRifDBfK4nILJo2w/3tJAunziDJ19HcXltbi3PnzgEAunbtCldXV7MF6yhVVVXw9PREZWUlz3cyA63eiEFvbMalOh2+fXAQhnT1Ex2JiNrAYJTx0JeHsSWjGH5ujvjp8aEI9bryOaPUfvVaA/q9vgl1WgNWPTaEE7CIyCzyyutww7+2QamQcOjFMfBxVYuOZHXaUhtc15mprq6uSEhIQEJCglUUTWR+G1ILcKlOhwB3Rwxkmx6R1VEqJLx/dxJigtxRWqPB3KWHUKPRi45lc7ZmFKNOa0CYtzM3CCciswn3cUHvEA8YjDI2nSoUHcfmcaQPtZssy/h8dxYAYGZyF6iU/HUiskZujip8Pqs//NwckVFYjaeXH4PByEl75vTz8XwApr2bJIltekRkPhPiggAAv6SxcOpo/KRL7XYwuxxp+VVwVCkwY1AX0XGI6DqEebvg05n9oFYpsDm9GP/cmCE6ks2oatBh22nTJNpbuektEZnZ+DjTROM9Z0tRWa8TnMa2sXCidvtidzYA4La+YeypJbIBfSO88fYdiQCAT3dmYcWh3BbuQa3x28kiaPVGdA9wQ0yQu+g4RGRjugW4oXuAG3QGmZvhdjAWTtQu58tqsanxH+fcYZFiwxCR2dyaGIKnRncHALy4Og37zpUJTmT9fm4c9c42PSLqKGzX6xwsnKhdluzJgSwDI3v6o1sAj6AS2ZKnx3THpMQQ6I0yHvn6CLJLa0VHslplNRrsOVsKgG16RNRxmtr1dmaWoJYDfjoMCydqs8p6HX44nAcAmNu4BwwR2Q5JkvDv2xPQJ9wLlfU6zF16CJV17Jtvjw1phTAYZSSEeSLSj9Nniahj9Ap2RxdfF2j0Rmw7XSw6js1i4URttuJQLmq1BvQMdMewbty3icgWOTko8enMfgjxdEJWaS0e/eYIdAaj6FhWZ23jpreTErjaREQdR5IkjGe7Xodj4URtojcYsWzveQDAnGGR7NcnsmEB7k744v4BcFUrsfdcGV766SSuY890u1NQWY+DOeUAgFsSgwWnISJbd3Nju962jGI06AyC09gmFk7UJhtPFiK/oh6+rmpM7hMqOg4RdbBewR54/64kSBLw3cFcLN6TIzqS1Vh3vAAAMDDSB8GezoLTEJGtSwjzRKiXM+q0BuzMLBEdxyaxcKI2aRpBfu/gLnByUApOQ0SdYUxsIF68uRcA4PX1p7A1g+NuW6N5ml4ftukRUceTJAnjepva9TayXa9DsHCiVjty/hKO5VZArVTg3sHc8JbInswdFoW7BoRDloEnvj2GjMIq0ZEsWnZpLVLzK6FUSLi58bwDIqKONiHe9HqzKd20fxyZFwsnarXFjatNk/uEwN/dUXAaIupMkiThtclxSI72Ra3WgLlLD6OkWiM6lsVa27jaNLSbH3zd+HpJRJ2jX4Q3/N0dUd2gx55zpaLj2BwWTtQqFy7V4Zc0U7/+3Bs4gpzIHqlVCiy8ty+i/FyRX1GPh786zEl7VyDLcnObHvduIqLOpFBIGNc7EACwMZXteubGwolaZdneHBhlYGg3X8QEeYiOQ0SCeLmo8cWs/vBwUuFobgU+3ZklOpLFySisxtniGqhVCtzU+AGGiKizTGicrvfbqULoeXDLrFg4UYtqNHosP8gNb4nIJNrfDa9O7g0AeH/zGWQWVQtOZFmaVptG9fSHh5OD4DREZG8GRfnA28UBl+p0OJhdLjqOTWHhRC364XAeqjV6RPu7YmSPANFxiMgCTOkTitExAdAajPjTD8d5VLORLMvN5zfdmsgtG4io86mUCoyNNa12czNc82LhRNdkMMpYvMc0FGLO0CgoFNzwlohMwyLeuC0eHk4qHL9Qic92ZYuOZBGO5VXgwqV6uKqVuDGGB5qISIwJ8aZ2vV9PFsJo5Mbl5sLCia5p06ki5JXXw8vFAdP6homOQ0QWJNDDCX+7JRYA8O7mTJwtZsvezymm1aaxsYFwVnOvOyISY2hXP7g7qVBcrcHR3Eui49gMFk50TU0jyGcMjOCHACL6g9v7hWFkT39o9Ub86ccTMNjxkU2DUcb6VNP00Vu56S0RCaRWKTCmF9v1zI2FE11V6oVKHMwph0ohYWZypOg4RGSBJEnCm7fFw91RhWO5Fc0HW+zRgawylFRr4OnsgGHd/EXHISI7N75x8+2NaYWQZfs9qGVOLJzoqr7YbRozfEtCMII8nQSnISJLFezpjL/e0gsA8PZvp3GupEZwIjGapundHB8EtYpvr0Qk1oge/nBRK5FfUY8TFypFx7EJfGWnKyqsbMC6E40b3g6LFpyGiCzdnf3DcUN3P2j0Rjxvhy17Wr2xuR1mEje9JSIL4OSgxKiepiE1bNczDxZOdEXL9uVAb5QxMMoH8WGeouMQkYWTJAlvTUuAm6MKR85fwtK9OaIjdapdZ0pQWa9DgLsjBkX5io5DRATg9+16BWzXMwMWTvQHdVo9vj2QC4Ab3hJR64V6OeMvN5ta9v79awZySmsFJ+o8TXs3TUwIhpLbNhCRhRgVEwC1SoGcsjpkFHLy6fVi4UR/sPJoPirrdYjwcWmeyEJE1Bp3DwzH0G6+aNCZWvbsYf+Qeq0Bv50qAgDcyjY9IrIgbo4qDO9uGlbDdr3rx8KJLmM0yljSOBVr9tBIHjklojaRJAlv3ZYAF7USB3PK8eW+HNGROtyvJwtRpzUg3McZfcK9RMchIrrMzfH/a9ej68PCiS6zPbMYWaW1cHdU4Y7+4aLjEJEVCvdxwYIJMQCAf248jdyyOsGJOtY3B84DAO7sFw5J4sEmIrIso3sFwkEpIbOoxm6nnpoLCye6zOe7TKtNdw+KgJujSnAaIrJW9wzqgsHRPqjXGfD8yuM227KXWVSNQzmXoFRIuHMADzYRkeXxdHbAkK5+AEx7OlH7sXCiZqcuVmHvuTIoFRJmDYkUHYeIrJhCIeFf0xLh7KDE/qxyfHMwV3SkDtE0SGdsr0AEenC/OyKyTBMap+ttSGW73vVg4UTNFu8xrTaNjwtCqJez4DREZO0ifF3w5/E9AQBvbkhHXrlttezVaw1YefQCAGDGoAjBaYiIrm5sbCAUEnDyYpXNt093JBZOBAAorm7AzymmcbocQU5E5jIzORIDI31QpzXghVUnbGofkbUnLqK6QY8IHxcM6+YnOg4R0VX5uv1vj7mNJ7nq1F4snAgA8PX+XGgNRiRFeKFvhLfoOERkIxQKCf+6PQFODgrsOVuG7w7miY5kNt80tundPTACCk4gJSILN6Fxuh7HkrcfCydCg86Ar/ebpkI9MCxacBoisjWRfq740zjTlL03NqQjv6JecKLrl5ZfieN5FXBQSrijf5joOERELRrX21Q4HcutQEGl9b8Oi8DCibDmWD7Ka7UI9XLGuN7c8JaIzO/+IZHo18UbNRo9Xlhp/S173zYOuxgfFww/N0fBaYiIWhbo4YR+XUxdRb9y1aldWDjZOVmWm4dC3D8kEiolfyWIyPyUjS17jioFdp0pxfeHrbdlr0ajx0/H8gEAMwZyKAQRWY+m6Xps12sffkq2c7vOlCKzqAauaiWmD+QeJETUcbr6u+HZm3oAAF5fl261rSI/peSjVmtAtL8rBkf7iI5DRNRq4xsLp0M55Sit0QhOY31YONm5L3abVpvu6B8ODycHwWmIyNbNHRaNpAgvVGv0WLAq1epa9mRZbt67acbACEgSh0IQkfUI83ZBQpgnjDLw28ki0XGsDgsnO7Ytoxg7MkugkIDZQyNFxyEiO6BUSPj37QlQqxTYfroEK4/mi47UJscvVOLkxSqoVQrc3o9DIYjI+oxvbtfjWPK2YuFkp6obdPjL6lQAwJyhUeji6yo4ERHZi24B7nhmjKll77W1J62qXeSbxgmkt8QHw8tFLTgNEVHbTYgLBgDsO1eGijqt4DTWhYWTnXrzlwwUVDagi68Lnr2pp+g4RGRnHrwhCnGhHqhq0OPNDRmi47RKZb0Oa0+YNgq/ZzCHQhCRdYryc0VMkDv0RhmbTrFdry1YONmhvedKm3v037otAc5qpeBERGRvVEoF/j45DpIErDx6AYdyykVHatHqoxfQoDOiZ6A7NwonIqvW1K63kdP12oSFk52p0+rxwkpTi949gyKQ3NVXcCIisldJEd64a4BpmudfV6dBZzAKTnR1siw37910z2AOhSAi69bUrrfrTCmqG3SC01gPFk525p3fMpFbXocQTye8MCFGdBwisnPPj4uBt4sDThdVY9neHNFxrurw+UvILKqBs4MSU5JCRcchIrouPQLdEO3nCq3BiK0ZxaLjWA0WTnbkyPlLzZvdvnFbPNw5fpyIBPN2VePP400Hcd7dlInCygbBia6sqb351sQQbt1ARFZPkiRMiGe7XluxcLITDToDnv/xOGQZmNY3DCN7BoiOREQEALizfziSIrxQqzXg9fWnRMf5g/JaLdanmsb2cigEEdmKpna97adLUK81CE5jHVg42Yn/bj2DcyW18HNzxN9u6SU6DhFRM4VCwt8nx0EhAetOFGDP2VLRkS6z8sgFaPVGxIV6ICHMS3QcIiKz6B3igTBvZ9TrDNiRyXa91mDhZAfS8iuxaEcWAOD1KXHce4SILE5cqCdmJkcCAP72Uxo0ess4+nnZUIhBXQSnISIyH0mSMKFxut6GVLbrtQYLJxunMxjxpx9PwGCUMTE+uHn8JBGRpZl/Uw/4uTkiq6QWn+/KFh0HgGmDyOzSWrg5qnBrYojoOEREZjW+sV1va0axxRywsmQsnGzcou3nkF5QBW8XB7xya2/RcYiIrsrDyQEvTjQNivjv1jO4cKlOcCLgm8bVpilJIXB1VAlOQ0RkXknhXgj0cESNRo/dZyyrTdoSsXCyYZlF1fjv1rMAgJcn9Ya/u6PgRERE1zalTygGRfmgQWfEq2vFDoooqdbg18ZpUzMGsk2PiGyPQiFhfG9TN9K6EwWC01g+Fk42ymCU8fyPJ6A1GDE6JgCT+7DFhIgsnyRJ+PuUOKgUEjadKsLWjCJhWX44kge9UUZShBdiQzyE5SAi6khNe9OtP1GA4mrL3BLCUrBwslFL9mQjJa8C7o4q/GNqPHe5JyKr0SPQHXOHRQEAXv75JBp0nd93bzTKzXs3cSgEEdmypAhvJEV4QWsw4uv9uaLjWDQWTjYop7QW//71NADgxYm9EOTpJDgREVHbPDm6O4I9nZBXXo+Pt5/r9OffeaYEFy7Vw8NJhVsSgjv9+YmIOlPTwapv9p8XcrDKWrBwsjFGo4w/rzwBjd6Iod18MX1AuOhIRERt5uqowt9uiQUALNpxDjmltZ36/E2rTdP6hcHJQdmpz01E1NnG9w5CqJczymq1+CklX3Qci8XCycZ8czAXB7LL4eygxFu3JbBFj4is1oS4INzQ3Q9avREv/3wSsix3yvMWVjZgS4ZpM8h7BkV0ynMSEYmkUiowa4ipLfmL3dmd9nprbVg42ZD8inq8tSEdAPD8+J4I93ERnIiIqP0kScJrk+OgViqwI7MEv57snA0aVxzKg8EoY2CUD7oFuHfKcxIRiTZ9QARc1EpkFtVg91mOJr8SFk42QpZlLFiVilqtAf27eGNWcqToSERE1y3KzxWPjIgGALy69hRqNfoOfT69wYjlh5qGQnC1iYjsh6ezA+7oFwbAtOpEf8TCyUasPJqPnZklUKsU+OftCVAo2KJHRLbhsVHdEO7jjILKBnyw9UyHPtf20yUoqGyAj6sa4+OCOvS5iIgszeyhUZAk02vh2eJq0XEsDgsnG1Bc1YDX1p4EADwzpge6+rsJTkREZD5ODkq8Mqk3AOCLXdk4U9Rxb+bfHDgPALi9XxgcVRwKQUT2JdLPFaNjAgEAi/fkiA1jgVg4WTlZlvHXNWmoatAjPtQTD94QJToSEZHZje4ViDG9AqE3yvjbT2kdcuJyXnkdtmeWAADuHsg2PSKyT02jyVcdvYBLtVrBaSwLCycrtz61AL+dKoJKIeFftydApeT/UiKyTS9PioWTgwL7s8rx8/GLZn/8FYfyIMvAsG5+iPJzNfvjExFZg8HRPogN9kCDzohvD3JD3N/jp2wrVl6rxcs/mVr0HhvVDb2CPQQnIiLqOOE+Lnh8VDcAwOvr01HVoDPbY+sMRqw4nAcAmMGhEERkxyRJal51WrY3B1q9UXAiy8HCyUrJsoyXfkpDWa0WPQPdmz9MEBHZsgeHRyPazxUl1Rq8uynTbI+7+VQRSqo18Hd3xNjYQLM9LhGRNZqUGAJ/d0cUV2uwIbVAdByLwcLJSv1nUybWnSiAQgL+dXsC1Cr+ryQi2+eoUuK1yXEATEdCT16sNMvjfnPA1I5yZ/8wOLDlmYjsnFqlwMzB3BD3/7OId4ePPvoIkZGRcHJywqBBg3Dw4MGr3nbp0qWQJOmyi5OTUyemFe+rfTn479azAIDXp8QjMdxLbCAiok40rLsfJiYEwygDf1uTBqPx+t7Qc0prsftsKSQJuGsA2/SIiADgnsFd4KhSIDW/EodyLomOYxGEF04rVqzA/Pnz8fLLL+Po0aNITEzEuHHjUFxcfNX7eHh4oKCgoPly/vz5Tkws1obUArz0s+m8pqfHdGcvPhHZpb9NjIWrWomjuRV4b3MmdmSW4GB2OVIvVOJscTUuXKpDWY0GtRo9DC0UVt81nvw8ooc/wn1cOiM+EZHF83FV47a+oQCAL3ZnCU5jGVSiA/znP//Bgw8+iNmzZwMAFi1ahPXr12Px4sV44YUXrngfSZIQFGR/GxPuO1eGp5enQJZNO9o/Nbq76EhEREIEeTrh6TE98I8N6figcQX+WtQqBZwdlKaLWglHlQLOatPXJy6Y2v3uGdSlo2MTEVmVOUOj8N3BPPx2qgi5ZXWI8LXvg0tCCyetVosjR45gwYIFzdcpFAqMGTMG+/btu+r9ampq0KVLFxiNRvTt2xdvvPEGevfufcXbajQaaDSa5q+rqqrM9wN0olMXq/DQl4ehNRgxvncQXpscB0mSRMciIhLm/qGRyCqtxenCKtTrjGjQGdCgM6BeZ0C91gDN7yZBafVGaPVGVNZfeRJfqJczRvX076zoRERWoXugO4b38MfOzBIs2ZuNlydd+fO2vRBaOJWWlsJgMCAw8PIJRoGBgcjIyLjifXr27InFixcjISEBlZWVePvttzFkyBCcPHkSYWFhf7j9m2++iVdffbVD8neWvPI6zFpyENUaPQZG+eC9u/pAqWDRRET2zUGpwJu3xV/1+0ajjAa9AQ06Y3Mx1VRYNTR+Xa8zQKMzYmCUD/fBIyK6grnDorAzswTfH8rDM2N7wMPJQXQkYYS36rVVcnIykpOTm78eMmQIevXqhU8++QR///vf/3D7BQsWYP78+c1fV1VVITw8vFOymkNZjQYzFx9ESbUGMUHu+Gxmfzg5KEXHIiKyeAqFBBe1Ci5q0UmIiKzX8O5+6B7ghjPFNfj+UB4euCFadCRhhB5e8/Pzg1KpRFFR0WXXFxUVtfocJgcHByQlJeHs2Sv3uDs6OsLDw+Oyi7Wo1egxZ+khZJfWItTLGcvmDISns/1W+URERETUuSRJwpzGDXGX7MmB3mC/G+IKLZzUajX69euHLVu2NF9nNBqxZcuWy1aVrsVgMCA1NRXBwcEdFVMIncGIR785iuMXKuHt4oBlcwYi0MO+xq4TERERkXhTk0Lh46pGfkU9fjtV1PIdbJTwhu758+fjs88+w7Jly5Ceno5HH30UtbW1zVP2Zs6cednwiNdeew2//fYbsrKycPToUdx77704f/48HnjgAVE/gtkZjTL+/OMJ7MwsgbODEovvH4BuAW6iYxERERGRHXJyUOKexi1wvtidLTiNOMLPcZo+fTpKSkrw0ksvobCwEH369MHGjRubB0bk5uZCofhffXfp0iU8+OCDKCwshLe3N/r164e9e/ciNjZW1I9gdv/cmIFVx/KhVEj4+N6+SIrwFh2JiIiIiOzYfYO7YNGOczhy/hJS8irQJ9xLdKROJ8myfH1brluZqqoqeHp6orKy0iLPd/p8VxZeX58OAHj7jkTc3u+PkwKJiIiIiDrb/BUpWHUsH5MSQ/Dfu5NExzGLttQGwlv16H/WHMtvLpr+PD6GRRMRERERWYymIREbUgtwsaJecJrOx8LJQuzMLMFzPxwHAMweGolHRtjvqEciIiIisjxxoZ4YFOUDg1HGsn05ouN0OhZOFuDEhQo88vUR6I0yJiWG4G8TYyFJ3OCWiIiIiCzL3MZVp+8O5KJOqxecpnOxcBIsu7QWs5ccQp3WgKHdfPH2HQlQKFg0EREREZHlGd0rEF18XVDVoMfKIxdEx+lULJwEKq5uwMzFB1BWq0VcqAcW3dsPjiql6FhERERERFekVEiYPSQSALB4Tw6MRvuZM8fCSZDqBh3uX3wIeeX16OLrgiX3D4S7k4PoWERERERE13RH/3C4O6mQXVqLbaeLRcfpNCycBKhu0GHussM4VVAFPzc1vpwzEP7ujqJjERERERG1yNVRhbsH2t+GuCycOllJtQZ3fbofB7PL4eaowtLZA9HF11V0LCIiIiKiVps1JBJKhYS958pw6mKV6DidgoVTJ8orr8Mdi/bi5MUq+LqqsfyhwYgL9RQdi4iIiIioTUK9nDE+LggAsHiPfaw6sXDqJBmFVZi2cC9yyuoQ5u2MHx8dwqKJiIiIiKxW02jyn1Muori6QXCajsfCqRMcyinHnYv2obhag56B7lj56BBE+bE9j4iIiIisV98IbyRFeEFrMOLr/bmi43Q4Fk4dbEt6Ee79/ACqGvTo38Ub3z+cjEAPJ9GxiIiIiIiuW9Oq0zf7z6NBZxCcpmOxcOpAK49cwENfHYFGb8SNMQH4au4geLpw5DgRERER2YbxvYMQ6uWMslotfrTxDXFZOHWQz3dl4dkfjsNglHFb31B8cl8/OKu5uS0RERER2Q6VUoHZQyMBAK+vP4VjuZfEBupALJzMTJZlvPVLBl5fnw4AeGBYFN6+PREOSv5VExEREZHtuX9IJEb19EeDzogHlh1Gblmd6Egdgp/mzUhvMOKFlalYtOMcAODP42Pw4sReUCgkwcmIiIiIiDqGSqnAhzP6oneIB8pqtbh/6UFU1GlFxzI7Fk5m0qAz4LFvjmLF4TwoJOCf0+Lx6MiukCQWTURERERk21wdVVh8/wAEezohq6QWD391BBq9bQ2LYOFkBlUNOsxafBC/nSqCWqXAwnv7YfqACNGxiIiIiIg6TaCHE5bMHgA3RxUOZJfjzz+egCzLomOZDQun61RSrcFdn+zHgexyuDmqsGz2QIzrHSQ6FhERERFRp4sJ8sDCe/tCpZCwJuUi3t2UKTqS2bBwug65ZXW4fdFenCqogp+bGssfGozkrr6iYxERERERCXNDd3+8MTUeAPDB1rP4/nCe4ETmwcKpndILqjBt0V6cL6tDuI8zfnxkCOJCPUXHIiIiIiIS7s4B4Xjixm4AgL+sSsXuM6WCE10/Fk5tJMsytmUU485P9qGkWoOYIHesfGQIIv1cRUcjIiIiIrIY88f2wOQ+IdAbZTz69RFkFFaJjnRdWDi1kt5gxNrjF3Hrh3swe+khVDfoMSDSGyseTkaAh5PoeEREREREFkWSJPzr9gQMjPJBtUaPOUsOoaiqQXSsdmPh1IJ6rQFf7svBqHe244nvjiE1vxJODgrcPyQSX84ZBE9nB9ERiYiIiIgskqNKiU/v64dof1dcrGzAnKWHUKvRi47VLpJsSzMCW6Gqqgqenp6orKyEh4fHVW9XXqvFl/tysGxvDi7V6QAA3i4OmDUkEjOTI+Hjqu6syEREREREVi23rA5TP96DslotbowJwKf39YNKKX4Np7W1AcDC6Q/fzyuvw+e7srDicB4adEYAQLiPMx68IRp39AuHs1rZ2ZGJiIiIiKzesdxLuOvT/dDojbh3cAT+PjkOkiQJzdSWwknVSZksXuqFSnyy8xw2pBbA2FhKxod64uER0RjfO8giKmIiIiIiImuVFOGN9+/qg0e/OYqv9+eii48rHhweLTpWq9l14STLMnadKcUnO89hz9my5uuH9/DHI8OjkdzVV3gVTERERERkK8bHBePFm3vh9fXp+MeGdIR6O+Pm+GDRsVrFbgun9Scu4ssjJUgvMI1FVCok3JoYggdviEZsyLWX6YiIiIiIqH3mDotCXnkdlu07j2dWpCDQwwn9uniLjtUiuz3HKfzp76FwdIGLWom7BkRgzrBIhHm7iI5HRERERGTzDEYZD391GJvTi+Hjqsbqx4agi2/n74valnOc7PbEHV9XBzx3Uw/sfeFGvDQplkUTEREREVEnUSokfHB3EuJDPVFeq8XsJYdwqVYrOtY12e2KU1FpOQJ8LX9JkIiIiIjIVhVXNWDqx3uRX1GPAZHe+HLOoE6dYs0Vp1ZwcuBYcSIiIiIikQI8nLBk9gC4O6pwKOcShry1Bf9YfwpZJTWio/2B3a44taaqJCIiIiKijrc/qwzzV6TgYmVD83VDuvpixqAI3BQbBLWqY9Z7uAHuNbBwIiIiIiKyPAajjO2ni/HNgVxsO12MpirFz02N2/uFY8bACET4mncuAQuna2DhRERERERk2fIr6rHiYC6WH8pDcbWm+fobuvvhnkERGN0rEA7K61+FYuF0Dc1/OSUXWTgREREREVkwncGI7adLsOJwHvacLW2+3t/NEdP6heGO/mEI8XRu9+NXVVXB0z+kVYWT3W6Ai3d6Ao6S6BRERERERHQVDgDGNl7g9Ltv6AEcaLxcD03r15DsdqoeERERERFRa9nvitOzpwG26hERERERWa3s0lr8cDgPq4/lo6JeBwC4o384Xru1d+seoKoKeCukVTe133OcOByCiIiIiMgmNOgMWHU0H39ZnQpJAn6aNxQJYV4t3o8b4BIRERERkd1wclBixqAITOkTAlkGXv75JIxG864PsXAiIiIiIiKbsODmXnBVK3EstwIrj14w62OzcCIiIiIiIpsQ6OGEJ0d3BwD8c2MGqhp0ZntsFk5ERERERGQzZg+NQrS/K0prtHhv0xmzPS4LJyIiIiIishlqlQKvTDJN1Vu2LweZRdVmeVwWTkREREREZFOG9/DHTbGBMBhlvPzTSZhjkDgLJyIiIiIisjl/uyUWjioF9mWVYUNq4XU/HgsnIiIiIiKyOeE+LnhkRFcAwOvrT6FOq7+ux2PhRERERERENunRkV0R5u2MgsoGfLzt3HU9FgsnIiIiIiKySU4OSvx1YiwA4NOdWcgprW33Y7FwIiIiIiIimzWudyBu6O4HrcGIv6871e7HYeFEREREREQ2S5IkvHJrbzgoJWzJKMbWjKJ2PQ4LJyIiIiIismld/d0wZ2gUAODVtafQoDO0+TFYOBERERERkc17YnR3BLg74nxZHb7Ynd3m+7NwIiIiIiIim+fmqMJfbu4FAPhw61lcrKhv0/1ZOBERERERkV2Y3CcEAyK9Ua8z4B8b0tt0XxZORERERERkFyRJwqu3xkEhAetPFODAubJW35eFExERERER2Y3YEA/cO7gLAODNja1fdWLhREREREREdmX+2B7wdnHA2eLWb4jLwomIiIiIiOyKl4saz4+PadN9WDgREREREZHdubN/OGKDPVp9exZORERERERkd5QKCX+Z2PpVJxZORERERERkl/qEe7f6thZROH300UeIjIyEk5MTBg0ahIMHD17z9j/88ANiYmLg5OSE+Ph4bNiwoZOSEhERERGRPRJeOK1YsQLz58/Hyy+/jKNHjyIxMRHjxo1DcXHxFW+/d+9e3H333Zg7dy6OHTuGKVOmYMqUKUhLS+vk5EREREREZC8kWZZlkQEGDRqEAQMG4MMPPwQAGI1GhIeH44knnsALL7zwh9tPnz4dtbW1WLduXfN1gwcPRp8+fbBo0aIWn6+qqgqenp6orKyEh0frTwYjIiIiIiLb0pbaQOiKk1arxZEjRzBmzJjm6xQKBcaMGYN9+/Zd8T779u277PYAMG7cuKveXqPRoKqq6rILERERERFRWwgtnEpLS2EwGBAYGHjZ9YGBgSgsLLzifQoLC9t0+zfffBOenp7Nl/DwcPOEJyIiIiIiuyH8HKeOtmDBAlRWVjZf8vLyREciIiIiIiIroxL55H5+flAqlSgqKrrs+qKiIgQFBV3xPkFBQW26vaOjIxwdHc0TmIiIiIiI7JLQFSe1Wo1+/fphy5YtzdcZjUZs2bIFycnJV7xPcnLyZbcHgE2bNl319kRERERERNdL6IoTAMyfPx+zZs1C//79MXDgQLz33nuora3F7NmzAQAzZ85EaGgo3nzzTQDAU089hREjRuCdd97BxIkTsXz5chw+fBiffvqpyB+DiIiIiIhsmPDCafr06SgpKcFLL72EwsJC9OnTBxs3bmweAJGbmwuF4n8LY0OGDMG3336Lv/71r/jLX/6C7t27Y82aNYiLixP1IxARERERkY0Tvo9TZ+M+TkREREREBFjRPk5ERERERETWgIUTERERERFRC1g4ERERERERtYCFExERERERUQtYOBEREREREbWAhRMREREREVELWDgRERERERG1gIUTERERERFRC1g4ERERERERtYCFExERERERUQtUogN0NlmWAQBVVVWCkxARERERkUhNNUFTjXAtdlc4lZWVAQDCw8MFJyEiIiIiIktQXV0NT0/Pa97G7gonHx8fAEBubm6LfzlEIlRVVSE8PBx5eXnw8PAQHYfoivh7SpaOv6Nk6fg7ahlkWUZ1dTVCQkJavK3dFU4Khem0Lk9PT/6SkkXz8PDg7yhZPP6ekqXj7yhZOv6OitfaxRQOhyAiIiIiImoBCyciIiIiIqIW2F3h5OjoiJdffhmOjo6ioxBdEX9HyRrw95QsHX9HydLxd9T6SHJrZu8RERERERHZMbtbcSIiIiIiImorFk5EREREREQtYOFERERERETUAhZORERERERELbC7wumjjz5CZGQknJycMGjQIBw8eFB0JKJmO3fuxKRJkxASEgJJkrBmzRrRkYiavfnmmxgwYADc3d0REBCAKVOm4PTp06JjEV1m4cKFSEhIaN5UNDk5Gb/88ovoWERX9NZbb0GSJDz99NOio1Ar2FXhtGLFCsyfPx8vv/wyjh49isTERIwbNw7FxcWioxEBAGpra5GYmIiPPvpIdBSiP9ixYwfmzZuH/fv3Y9OmTdDpdLjppptQW1srOhpRs7CwMLz11ls4cuQIDh8+jBtvvBGTJ0/GyZMnRUcjusyhQ4fwySefICEhQXQUaiW7Gkc+aNAgDBgwAB9++CEAwGg0Ijw8HE888QReeOEFwemILidJElavXo0pU6aIjkJ0RSUlJQgICMCOHTswfPhw0XGIrsrHxwf//ve/MXfuXNFRiAAANTU16Nu3Lz7++GO8/vrr6NOnD9577z3RsagFdrPipNVqceTIEYwZM6b5OoVCgTFjxmDfvn0CkxERWafKykoApg+lRJbIYDBg+fLlqK2tRXJysug4RM3mzZuHiRMnXva5lCyfSnSAzlJaWgqDwYDAwMDLrg8MDERGRoagVERE1sloNOLpp5/G0KFDERcXJzoO0WVSU1ORnJyMhoYGuLm5YfXq1YiNjRUdiwgAsHz5chw9ehSHDh0SHYXayG4KJyIiMp958+YhLS0Nu3fvFh2F6A969uyJlJQUVFZW4scff8SsWbOwY8cOFk8kXF5eHp566ils2rQJTk5OouNQG9lN4eTn5welUomioqLLri8qKkJQUJCgVERE1ufxxx/HunXrsHPnToSFhYmOQ/QHarUa3bp1AwD069cPhw4dwvvvv49PPvlEcDKyd0eOHEFxcTH69u3bfJ3BYMDOnTvx4YcfQqPRQKlUCkxI12I35zip1Wr069cPW7Zsab7OaDRiy5Yt7HsmImoFWZbx+OOPY/Xq1di6dSuioqJERyJqFaPRCI1GIzoGEUaPHo3U1FSkpKQ0X/r374977rkHKSkpLJosnN2sOAHA/PnzMWvWLPTv3x8DBw7Ee++9h9raWsyePVt0NCIApik7Z8+ebf46OzsbKSkp8PHxQUREhMBkRKb2vG+//RY//fQT3N3dUVhYCADw9PSEs7Oz4HREJgsWLMCECRMQERGB6upqfPvtt9i+fTt+/fVX0dGI4O7u/ofzQl1dXeHr68vzRa2AXRVO06dPR0lJCV566SUUFhaiT58+2Lhx4x8GRhCJcvjwYYwaNar56/nz5wMAZs2ahaVLlwpKRWSycOFCAMDIkSMvu37JkiW4//77Oz8Q0RUUFxdj5syZKCgogKenJxISEvDrr79i7NixoqMRkZWzq32ciIiIiIiI2sNuznEiIiIiIiJqLxZORERERERELWDhRERERERE1AIWTkRERERERC1g4URERERERNQCFk5EREREREQtYOFERERERETUAhZORETUKbZv3w5JklBRUdHpzy1JEiRJgpeX13U/Vmt+jqVLl173c40cObI5d0pKynU9FhERXT8WTkREZHYjR47E008/fdl1Q4YMQUFBATw9PYVkWrJkCTIzMzvluaZPn37Zc73yyivo06dPmx5j1apVOHjwoJmTERFRe6lEByAiIvugVqsRFBQk7Pm9vLwQEBBwXY+h0+ladTtnZ2c4Oztf13P5+Pigqqrquh6DiIjMhytORERkVvfffz927NiB999/v7nVLCcn5w8tbk3tbOvWrUPPnj3h4uKC22+/HXV1dVi2bBkiIyPh7e2NJ598EgaDofnxNRoNnnvuOYSGhsLV1RWDBg3C9u3b25V14cKF6Nq1K9RqNXr27Imvvvrqsu9LkoSFCxfi1ltvhaurK/7xj380f2/Pnj1ISEiAk5MTBg8ejLS0tObv/b5Vb+nSpXj11Vdx/Pjx5r+PpUuXQpZlvPLKK4iIiICjoyNCQkLw5JNPtuvnICKijscVJyIiMqv3338fmZmZiIuLw2uvvQYA8Pf3R05Ozh9uW1dXhw8++ADLly9HdXU1brvtNkydOhVeXl7YsGEDsrKyMG3aNAwdOhTTp08HADz++OM4deoUli9fjpCQEKxevRrjx49Hamoqunfv3uqcq1evxlNPPYX33nsPY8aMwbp16zB79myEhYVh1KhRzbd75ZVX8NZbb+G9996DSqVCVlYWAOBPf/oT3n//fQQFBeEvf/kLJk2ahMzMTDg4OFz2PNOnT0daWho2btyIzZs3AwA8PT2xcuVKvPvuu1i+fDl69+6NwsJCHD9+vE1/10RE1HlYOBERkVl5enpCrVbDxcWlxdY8nU7XvOoDALfffju++uorFBUVwc3NDbGxsRg1ahS2bduG6dOnIzc3F0uWLEFubi5CQkIAAM899xw2btyIJUuW4I033mh1zrfffhv3338/HnvsMQDA/PnzsX//frz99tuXFU4zZszA7Nmzm79uKpxefvlljB07FgCwbNkyhIWFYfXq1bjzzjsvex5nZ2e4ublBpVJd9veRm5uLoKAgjBkzBg4ODoiIiMDAgQNbnZ+IiDoXW/WIiEgYFxeX5qIJAAIDAxEZGQk3N7fLrisuLgYApKamwmAwoEePHnBzc2u+7NixA+fOnWvTc6enp2Po0KGXXTd06FCkp6dfdl3//v2veP/k5OTmP/v4+KBnz55/uO+13HHHHaivr0d0dDQefPBBrF69Gnq9vg0/ARERdSauOBERkTD/v61NkqQrXmc0GgEANTU1UCqVOHLkCJRK5WW3+32xZU6urq4d8rjh4eE4ffo0Nm/ejE2bNuGxxx7Dv//9b+zYseMPfwdERCQeV5yIiMjs1Gr1ZQMdzCUpKQkGgwHFxcXo1q3bZZe2Tuzr1asX9uzZc9l1e/bsQWxsbKvuv3///uY/X7p0CZmZmejVq9cVb3u1vw9nZ2dMmjQJH3zwAbZv3459+/YhNTW1DT8FERF1Fq44ERGR2UVGRuLAgQPIycmBm5sbfHx8zPK4PXr0wD333IOZM2finXfeQVJSEkpKSrBlyxYkJCRg4sSJrX6sP/3pT7jzzjuRlJSEMWPGYO3atVi1alXzAIeWvPbaa/D19UVgYCBefPFF+Pn5YcqUKVe8bWRkJLKzs5GSkoKwsDC4u7vju+++g8FgwKBBg+Di4oKvv/4azs7O6NKlS6t/BiIi6jxccSIiIrN77rnnoFQqERsbC39/f+Tm5prtsZcsWYKZM2fi2WefRc+ePTFlyhQcOnQIERERbXqcKVOm4P3338fbb7+N3r1745NPPsGSJUswcuTIVt3/rbfewlNPPYV+/fqhsLAQa9euhVqtvuJtp02bhvHjx2PUqFHw9/fHd999By8vL3z22WcYOnQoEhISsHnzZqxduxa+vr5t+jmIiKhzSLIsy6JDEBERdSRJkrB69eqrrghZqpycHERFReHYsWPo06eP6DhERHaNK05ERGQX7r77boSFhYmO0WoTJkxA7969RccgIqJGXHEiIiKbd/bsWQCAUqlEVFSU4DStk5+fj/r6egBARETEVdsAiYioc7BwIiIiIiIiagFb9YiIiIiIiFrAwomIiIiIiKgFLJyIiIiIiIhawMKJiIiIiIioBSyciIiIiIiIWsDCiYiIiIiIqAUsnIiIiIiIiFrAwomIiIiIiKgFLJyIiIiIiIha8H+/QDSADQAPxwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "ax = plt.subplot(111)\n",
    "ax.set_xlabel(\"time [orbits]\")\n",
    "ax.set_xlim([0,sim.t/(2.*np.pi)])\n",
    "ax.set_ylabel(\"distance\")\n",
    "plt.plot(times/(2.*np.pi), distances);\n",
    "plt.plot([0.0,12],[0.2,0.2]); # Plot our close encounter criteria;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We did indeed find the close encounter correctly. We can now search for the two particles that collided and, for this example, merge them. To do that we'll first calculate our new merged planet coordinates, then remove the two particles that collided from REBOUND and finally add the new particle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of particles at the beginning of the simulation: 3.\n",
      "Two particles had a close encounter (d<exit_min_distance).\n",
      "Number of particles at the end of the simulation: 2.\n"
     ]
    }
   ],
   "source": [
    "from itertools import combinations\n",
    "def mergeParticles(sim):\n",
    "    # Find two closest particles\n",
    "    min_d2 = 1e9 # large number\n",
    "    ps = sim.particles\n",
    "    for i1, i2 in combinations(range(sim.N),2): # get all pairs of indices\n",
    "        dp = ps[i1] - ps[i2]   # Calculates the coponentwise difference between particles \n",
    "        d2 = dp.x*dp.x+dp.y*dp.y+dp.z*dp.z\n",
    "        if d2<min_d2:\n",
    "            min_d2 = d2\n",
    "            col_i1 = i1\n",
    "            col_i2 = i2\n",
    "    \n",
    "    cp1 = ps[col_i1]\n",
    "    cp2 = ps[col_i2]\n",
    "    # Merge two closest particles\n",
    "    \n",
    "    sum_mass = cp1.m + cp2.m\n",
    "    mergedPlanet = (cp1*cp1.m + cp2*cp2.m)/sum_mass \n",
    "    mergedPlanet.m  = sum_mass\n",
    "    sim.remove(index=col_i2) # Note: Removing a particle changes the sim.particles\n",
    "    sim.remove(index=col_i1) #       array and the particle indicies.\n",
    "    sim.add(mergedPlanet, assignHash=True)\n",
    "\n",
    "sim = setupSimulation() # Resets everything\n",
    "sim.exit_min_distance = 0.15\n",
    "print(\"Number of particles at the beginning of the simulation: %d.\"%sim.N)\n",
    "for i,time in enumerate(times):\n",
    "    try:\n",
    "        sim.integrate(time)\n",
    "    except rebound.Encounter as error:\n",
    "        print(error)\n",
    "        mergeParticles(sim)\n",
    "print(\"Number of particles at the end of the simulation: %d.\"%sim.N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "We can achieve the same outcome by using more of the built-in functionality of REBOUND. For that, we set the radius of the particles to their Hill radius. In practice, you might want to use the physical radius, but for this example, we want the collision to occur in a short amount of time and therefore inflate the particle radii. We set the collision detection routine to `direct` which will do a $O(N^2)$ collision search between all particles. The `collisions_resolve` call-back function is set to `merge`, which will merge the particles together, assuming mass and momentum conservation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setupSimulation():\n",
    "    sim = rebound.Simulation()\n",
    "    sim.integrator = \"ias15\" # IAS15 is the default integrator, so we don't really need this line\n",
    "    sim.add(m=1.)\n",
    "    sim.add(m=1e-3, a=1., r=np.sqrt(1e-3/3.)) # we now set collision radii!\n",
    "    sim.add(m=5e-3, a=1.25, r=1.25*np.sqrt(5e-3/3.))\n",
    "    sim.move_to_com()\n",
    "    return sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Particles in the simulation at t=   0.0: 3\n",
      "Particles in the simulation at t= 100.0: 2\n"
     ]
    }
   ],
   "source": [
    "sim = setupSimulation()\n",
    "sim.collision = \"direct\"\n",
    "sim.collision_resolve = \"merge\"\n",
    "\n",
    "print(\"Particles in the simulation at t=%6.1f: %d\"%(sim.t,sim.N))\n",
    "sim.integrate(100.)\n",
    "print(\"Particles in the simulation at t=%6.1f: %d\"%(sim.t,sim.N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the built-in collision detection and apply our own function to resolve the collision. By default, if we don't set the sim.collision function pointer, `REBOUND` will raise a `Collision` exception when a collision occurs, which we can catch. \n",
    "\n",
    "An indirect way of checking which particles collided is to check which ones have a `last_collision` time equal to the current simulation time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Particles [1, 2] collided\n"
     ]
    }
   ],
   "source": [
    "sim = setupSimulation()\n",
    "sim.collision = \"direct\"\n",
    "# we don't set sim.collision_resolve this time\n",
    "\n",
    "try:\n",
    "    sim.integrate(100.)\n",
    "except rebound.Collision:\n",
    "    collided = []\n",
    "    for p in sim.particles:\n",
    "        if p.last_collision == sim.t:\n",
    "            collided.append(p.index)\n",
    "    # Custom resolution\n",
    "\n",
    "print(\"Particles {0} collided\".format(collided))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
